home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / basic / ubmalm.zip / dldsys.ub < prev    next >
Text File  |  1990-08-22  |  5KB  |  127 lines
  1.    10   ' Driver program for LDIOSYS
  2.    11   ' Gets input from a file.
  3.    12   ' 19 June 1990
  4.    15   word 20
  5.    20   dim A(15,16),V(15,15)
  6.    22   input ="Bexample.ubd"
  7.    24   input R:input C
  8.    26   for I=1 to R:for J=1 to C+1:input A(I,J):next J:next I
  9.    28   input =input
  10.   100   gosub *Ldiosys(R,C,&A(),&V(),&Soln,&Rank)
  11.   290   print =print + lprint
  12.   300   gosub *Report(R,C,&A(),&V(),&Soln,&Rank)
  13.   302   print:print:print:print =lprint +"Malm.UBD"
  14.   305   gosub *Report2(R,C,&A(),&V(),&Soln,&Rank)
  15.   306   print:print:print:print:print:
  16.   310   print =print
  17.   400   end
  18.   490   '
  19.   495   '
  20.   800   *Report(R,C,&A(),&V(),&Soln,&Rank)
  21.   810   ' Outputs the results in a format that is O.K. if the matrices are not
  22.   820   ' too big nor are the numbers.
  23.   830   ' 18 June 1990
  24.   850   local I%,J%,S
  25.   860   if Soln=0 then print "No solution." endif
  26.   870   if Rank=0 then print "The rank is zero." else
  27.   875   :print "The rank is ";Rank
  28.   880   :print "The augmented coeff. matrix followed by the record matrix:"
  29.   890   :for I%=1 to R:print "        ";
  30.   900   :for J%=1 to C+1:print A(I%,J%),
  31.   910   :next J%:print:next I%
  32.   920   :for I%=1 to C:print "        ";
  33.   930   :for J%=1 to C:print V(I%,J%),
  34.   940   :next J%:print:next I%
  35.   950   :for I%=1 to 3:print:next I%
  36.   960   :if Soln then
  37.   970   ::for I%=1 to C:S=0:print "        ";
  38.   980   ::for J%=1 to Rank:S+=V(I%,J%)*A(J%,C+1):next J%
  39.   990   ::print S,:print "          ";
  40.  1000   ::for J%=Rank+1 to C:print V(I%,J%),:next J%
  41.  1010   ::print:next I%
  42.  1020   ::endif ' if soln
  43.  1030   :endif:return ' End of SubroutineReport.
  44.  1090   '
  45.  1095   '
  46.  1100   *Report2(R,C,&A(),&V(),&Soln,&Rank)
  47.  1110   ' Outputs the results one number to a line.
  48.  1130   ' 19 June 1990.
  49.  1140   local I%,J%,S
  50.  1160   print "The augmented coefficient matrix by columns:":print
  51.  1170   for J%=1 to C+1:for I%=1 to R:print A(I%,J%):next I%:print:next J%
  52.  1200   print:print "The record matrix by columns:":print
  53.  1210   for J%=1 to C:for I%=1 to C:print V(I%,J%):next I%:print:next J%
  54.  1230   print:print "The rank is ";Rank
  55.  1240   if Soln=0 then print "No solution":return endif
  56.  1250   print "A particular solution (column):":print
  57.  1260   for I%=1 to C:S=0
  58.  1270   for J%=1 to Rank:S+=V(I%,J%)*A(J%,C+1):next J%
  59.  1280   print S:next I%:print
  60.  1300   print "The coeff. matrix of the parameters(by columns):":print
  61.  1310   for I%=Rank+1 to C:for J%=1 to C:print V(J%,I%):next J%:print:next I%
  62.  1320   return ' End of subroutine Report2.
  63.  1330   '
  64.  1335   '
  65.  2320   *LDiosys(R,C,&A(),&V(),&Soln,&Rank)
  66.  2330   ' Solve the linear Diophantine system represented by the (augmented)
  67.  2340   ' matrix A().  Matrix V() contains a record of the manipulations.
  68.  2350   ' A is R by C+1, while V is C by C.  Soln is 1 if there is a solution
  69.  2360   ' and is 0 if there is no solution.  Rank is the rank of the system.
  70.  2370   '
  71.  2380   ' NOTE - there is a GLOBAL variable J% - required to receive the value
  72.  2390   ' of the function Pivotind, which LDiosys calls.
  73.  2400   '
  74.  2410   ' 19 June 1990
  75.  2420   local K%,M%,N%,I%,Done,Te,Found,P
  76.  2430   for I%=1 to C:for K%=1 to C:V(I%,K%)=0:next K%:next I%
  77.  2440   for I%=1 to C:V(I%,I%)=1:next I%
  78.  2450   Done=0:I%=0:Soln=1
  79.  2460   while and{I%<R,I%<C,Done=0}
  80.  2470   inc I%:J%=fnPivotind(I%)
  81.  2480   if J%=0 then K%=I%:Found=0
  82.  2490   :while and{Found=0,K%<R}
  83.  2500   :inc K%:J%=fnPivotind(K%)
  84.  2510   :if J%<>0 then Found=1 endif
  85.  2520   :wend
  86.  2530   :if Found=0 then Done=1 else
  87.  2540   :for M%=1 to C+1:swap A(K%,M%),A(I%,M):next M%
  88.  2550   :endif endif
  89.  2560   if Done=0 then
  90.  2570   :repeat K%=J%:P=A(I%,J%)
  91.  2580   :for N%=I% to C
  92.  2590   :if N%<>J% then Te=A(I%,N%)\P
  93.  2600   :for M%=I% to R:A(M%,N%)-=Te*A(M%,J%):next M%
  94.  2610   :for M%=1 to C:V(M%,N%)-=Te*V(M%,J%):next M%
  95.  2620   :endif:next N%
  96.  2630   :J%=fnPivotind(I%)
  97.  2640   :until J%=K%
  98.  2650   :if I%<>J% then
  99.  2660   :for M%=I% to R:swap A(M%,I%),A(M%,J%):next M%
  100.  2670   :for M%=1 to C:swap V(M%,I%),V(M%,J%):next M%
  101.  2680   :endif
  102.  2690   :if A(I%,C+1)@A(I%,I%)<>0 then Done=1 else
  103.  2700   :A(I%,C+1)\=A(I%,I%):A(I%,I%)=1
  104.  2710   :for M%=I%+1 to R:A(M%,C+1)-=A(M%,I%)*A(I%,C+1):A(M%,I%)=0:next M%
  105.  2720   :endif
  106.  2730   :endif ' done=0
  107.  2740   wend
  108.  2750   Te=min(R,C):Rank=0
  109.  2760   for K%=1 to Te:if A(K%,K%)<>0 then inc Rank endif next K%
  110.  2770   K%=0
  111.  2780   while and{K%<Rank,Soln}:inc K%:if A(K%,K%)<>1 then Soln=0 endif wend
  112.  2790   for K%=Rank+1 to R:if A(K%,C+1)<>0 then Soln=0 endif:next K%
  113.  2800   return ' End of subroutine LDiosys.
  114.  2810   '
  115.  2820   fnPivotind(Kk%)
  116.  2830   ' Returns the column index of the smallest (in absolute value )
  117.  2840   ' non-zero entry in row Kk% of the coefficient matrix A.
  118.  2850   ' If all entries are zero, then Pivotind = 0.
  119.  2855   ' 15 June 1990
  120.  2860   local T,S,K%,J%
  121.  2870   T=abs(A(Kk%,1))
  122.  2880   if T<>0 then K%=1 else K%=0 endif
  123.  2890   for J%=2 to C
  124.  2900   S=abs(A(Kk%,J%))
  125.  2910   if and{S<>0,or{T=0,S<T}} then T=S:K%=J% endif
  126.  2920   next J%:return(K%) ' End of function Pivotind.
  127.